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Abstract We carry out numerical and mathematical investigations of shear Alfven waves inside of a spheri- 
cal shell filled with an incompressible conducting fluid, and bathed in a strong dipolar magnetic field. We fo- 
cus on axisymmetric toroidal and non-axisymmetric modes, in continuation of a previous work by Rincon & 
Rieutord (2003, A&A, 398, 663). Analytical expressions are obtained for toroidal eigenmodes and their corre- 
sponding frequencies at low diffusivities. These oscillations behave like magnetic shear layers, in which the mag- 
netic poles play a key role, and hence become singular when diffusivities vanish. It is also demonstrated that 
non-axisymmetric modes are split into two categories, namely poloidal or toroidal types, following similar asymp- 
totic behaviours as their axisymmetric counterparts when the diffusivities become arbitrarily small. 
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1. Introduction 

Numerous astrophysical systems exhibit a pulsating be- 
haviour that can be significantly influenced by the Lorentz 
force when a strong magnetic field is present. This may for 
instance be the case in neutron stars and magnetic white 
dwarfs (Lou, 1995). Planetary cores, which are known to 
sustain a strong dynamo (Stevenson, 1983), are also likely 
to fall into this category. 

One of the most exciting examples of couplings be- 
tween pulsation and magnetism is given by the seismo- 
logical activity of roAp stars. This class of stars, discov- 
ered by Kurtz (1978), exhibits several kG (almost) dipolar 
magnetic fields and luminosity variations on periods rang- 
ing from 5 to 15 min. These oscillations seem to be well 
approximated by a single spherical harmonic I = 1 lined 
up with the magnetic axis, suggesting a strong mixing be- 
tween high order p- modes and Alfvenic type oscillations. 

Many different models have been developed to obtain 
a satisfying picture of the asteroseismology of these stars. 
Following theoretical work by Biront et al. (1982), Roberts 
& Soward (1983) and Campbell & Papaloizou (1986), 
Dziembowski & Goode (1996) have studied acoustic star 
models enveloped by a layer in which magnetic effects be- 
come dominant. Using a boundary layer approximation, 
they came up with an outer boundary condition which 
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was then applied for the calculation of adiabatic acous- 
tic modes. This model has undergone a lot of refinements 
in order to take into account new physical processes. The 
latest version, proposed by Bigot & Dziembowski (2002), 
incorporates the centrifugal force (a non-axisymmetric ef- 
fect, since the rotational and magnetic axis are often tilted 
in roAp stars) and suggests that the axis of the modes is 
not necessarily lined up with the magnetic axis. In spite of 
these improvements, there are still non-negligible discrep- 
ancies between the magnetically shifted eigenfrequencies 
computed from these models and the observed ones, show- 
ing that the precise coupling mechanism occurring in the 
surface layers is likely to be more complex. An important 
point is that a single t value is sometimes assumed to 
be sufficient to describe the oscillations. This may not be 
the case, owing to the dipolar structure of the permanent 
magnetic field which induces a coupling between spherical 
harmonics thus producing a whole spectrum of ts (e.g. 
Rincon & Rieutord 2003). 

Motivated by the observation that chemical peculiar- 
ities are observed near the magnetic poles of roAp stars, 
Balmforth et al. (2001) have tried to determine what pre- 
cise physical phenomena were occurring in the polar and 
equatorial regions. Since the magnetic field is almost hor- 
izontal near the equator and vertical near the poles, con- 
vection is certainly inhibited in the latter region, allowing 
the diffusion of different chemical elements (an excess of 
helium is observed on the polar caps). Oscillations trig- 
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gered by a K-mechanism may therefore preferentially be 
observed in this region. This approach stresses the impor- 
tance of a global description of the eigenmodes in such 
stars. 

The work presented here originates partly from the 
preceding remarks. It continues the study described in 
Rincon & Rieutord 2003 (hereafter referred to as Paper I) 
and should be viewed as a preliminary step to obtain 
global models of the more complex magneto-gravito- 
acoustic oscillations. We aim at understanding some of 
the underlying physical mechanisms potentially involved 
in roAp stars or similar astrophysical objects, by using 
an approach that rigorously treats the couplings induced 
by the geometry of the magnetic field. Our highly sim- 
plified model consists of a non-rotating spherical shell of 
incompressible magnetised fluid bathed in a dipolar mag- 
netic field, with small magnetic and kinematic diffusivi- 
ties. As a consequence of incompressibility, the modes we 
compute are shear Alfven waves. In Paper I, a spherical 
harmonic decomposition of the linearised MHD equations 
was obtained and results regarding poloidal axisymmetric 
modes were presented (a short reminder of the classifi- 
cation of modes is given later on in the paper). It was 
shown that the least-damped modes are near the mag- 
netic poles and exhibit internal shear layers which can 
potentially play a role in mode selection. In the present 
study, we focus on the other types of shear Alfvenic oscil- 
lations in spherical shells, namely axisymmetric toroidal 
and non-axisymmctric modes, and characterise numeri- 
cally and mathematically their geometry, eigenspectrum, 
and some by-products like boundary layers. 

In Sect. 2, we recall the basic physical ingredients of 
the model and shortly describe our numerical strategy. 
Sect. 3 is devoted to the phenomenological description and 
analytical study of axisymmetric toroidal modes. Sect. 4 
covers non-axisymmetric modes and their resemblance to 
their axisymmetric counterparts. Finally, Sect. 5 and 6 
conclude the paper. Note that most of the mathematical 
details involved in Sect. 3 are given in App. A and App. B 
to preserve the clarity of the manuscript. 



2. Description of the model 

2.1. Basic equations 

We first give a brief description of the model that is used 
in calculating toroidal and non-axisymmetric modes and 
which was already used in Paper I. More details on the ba- 
sic equations and their expansion onto the spherical har- 
monic base are given in Paper I. 

The "star" we work with is a spherical shell of incom- 
pressible plasma of density p Ql with a radius R and an 
aspect ratio 77. In all the numerical examples and figures, 
we use 77 = 0.35 1 . Within this shell is a dipolar magnetic 
field generated in a perfectly conducting core: 



This is the aspect ratio of the Earth's liquid core or that 
of the radiative zone of a 3 M star. 
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We use the following underlined dimensionless variables: 



Rr, v = V A v, B = B B, 



R 



where Va, the Alfven velocity, is given by the following 
expression, \i D being the magnetic permeability: 

B n 



V A = 



y/PofJ-o 



(2) 



We apply the linearised magneto-hydrodynamic equations 
to the kinetic v and magnetic b perturbations (we drop 
the underlined notation). We assume these perturbations 
have a time dependence of the form e At , where A = r + iiv 
(r is the damping rate, uj the pulsation and i 2 = —1). This 
leads to the following set of equations: 



V-v =0, 
V • b =0, 

AV x v = V x ((V x b) x B) + EV x Av, 
Xb = V x (v x B) + E m Ab. 



(3) 



where B denotes the permanent dipolar magnetic field. 
The parameters E and E m are non-dimensional forms of 
the kinematic and magnetic diffusivities, respectively, and 
are given by the following expressions: 

(4) 



E 



RV A ' 



E n 



<J fJ, RV A 



where <r is the conductivity of the fluid. In the case of 
roAp stars, these parameters take on the following typical 
values : 



E„ 



10" 



£~ 10 



-13 



2.2. Boundary conditions 

The boundary conditions used in our model are as follows. 
For the velocity, we use stress-free, non-penetrative, con- 
ditions to minimise the importance of boundary layers. 
This leads to the following conditions, valid for r = 77 and 
r = 1: 



0, 



0. (5) 



or v r / or \ r 

As for the magnetic field, different conditions apply for the 
inner and outer boundaries. On the interior, the perturba- 
tion to the electric field is perpendicular to the conducting 
core, and the perturbation to the magnetic field is tangent. 
This gives the following three equations: 



0, 



Em d 
r or 

r or 



-v B r 



-v v B r . 



(6) 



When the other boundary conditions are taken into ac- 
count, and when A is different from zero, these equations 
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are not independent and only correspond to two condi- 
tions. On the outer boundary and beyond, the magnetic 
field is continuous and potential, since there are no cur- 
rents in empty space. By using the continuity of this field, 
and the idea that it does not diverge toward infinity, it 
is possible to define two boundary conditions, which are 
more easily expressed in spherical harmonics (see Eq. (8)). 

Eq. (3), together with these boundary conditions, de- 
fines a generalised eigenvalue problem, where A is the 
eigenvalue and (v, b) is the eigenvector which can be com- 
puted numerically. 

2.3. Harmonic projection 

In order to solve Eq. (3), the fields v and b are projected 
onto the harmonic base: 

oo I 

v = E E uiRT + viST + wiTf, 

e=om=-e rrj\ 

b = E E alRT + blST + ciTT, 

t=0 m=-t 

in which R™, S™, and T™ are the normalised spherical 
harmonics: 

Rf = Y e m e r , Sf = rVY e m , If = rV x R™. 

The full harmonic decomposition of Eq. (3) is given in 
App. A of Paper I. The outer boundary condition on the 
magnetic field reads: 

dot (e + 2)aj n _ 
dr + r ' (8) 



2.4. Classification and symmetries 

The various eigenmodes fit naturally into different cate- 
gories. Firstly, as was already shown in Paper I, there is no 
coupling between different m's. Hence, eigenmodes will be 
made up of only one to. This leads to two types of modes: 
axisymmetric oscillations (to = 0) and non-axisymmetric 
ones (to ytz 0) . Secondly, within the axisymmetric category, 
it is possible to distinguish between poloidal modes and 
toroidal ones, as the corresponding equations fully decou- 
ple for to = (see App. A of Paper I). Poloidal modes are 
made up of Uq, Vq, a^, and 6g functions which correspond 
to e r and eg components. Toroidal modes are made up of 
Wq and Cq functions and are in the e v direction. 

A certain number of symmetries are present in the 
physical system and lead to a few simplifications. As was 
pointed out in Paper I, a parity can be defined for eigen- 
modes. However, there was a slight confusion as to the 
parity of toroidal components (see Eq. (19) of Paper I) 
since a mode is called even when the velocity perturba- 
tion is symmetric with respect to the equator and the 
magnetic perturbation antisymmetric, and a mode is odd 



in the reverse situation. The corrected form for toroidal 
eigenvectors reads: 

(v,b) = {w 2k +\c 2k ) even, 

(v,b) = {wl k ,cf +1 ) odd. W 

Non-axisymmetric modes also have a global parity. In 
fact it is possible to anticipate this result without the use 
of spherical harmonics. If we define M. as being the op- 
erator that gives the mirror image of a vector field with 
respect to the equator 2 , and if (v, b, A) is a solution, then 
(A4v, —Mb, A) is also a solution. By combining these two 
solutions, it is possible to extract an even part and an odd 
part from the original solution. Eigenvectors will take on 
the following form: 

{v,b) = (u™ +2k ,w™ +2k+1 ,a™+ 2k+1 ,c™ +2fe ) even, 
(v,b) = «+ 2fe + 1 ,<+ 2fe ,a™+ 2fe ,C+ 2fe + 1 ) odd. 

(10) 

This symmetry enables us to work with half as many com- 
ponents for a given resolution, which is advantageous from 
the numerical point of view. Useful information about the 
eigenspectra can be deduced from the symmetry of B with 
respect to any meridian. For a given solution (y, 6, A) it is 
possible to create a second solution (Sv, Sb, A), where S is 
the operator 3 that gives the mirror image of a vector field 
with respect to the meridian that passes through <p = 0. 
When applied to spherical harmonics, the azimuthal or- 
der to changes to — m. From this we conclude that the 
eigenspectrum of modes with azimuthal order to is the 
same as that of —to. This also leads to the decoupling 
of poloidal and toroidal components in the axisymmetric 
case (the poloidal modes being symmetric with respect 
to S and the toroidal modes antisymmetric). By taking 
into account that B is real, we deduce that the eigenspec- 
trum for to is the conjugate of that for —to. Therefore, 
each eigenspectrum is symmetrical with respect to the real 
axis. Practically, this means that we only need to explore 
eigenvalues for positive azimuthal orders, and only on the 
upper half of the complex plane (ui > 0). 

If rotation were taken into account, a number of these 
symmetries would break down. Apart from the situation 
where the rotation and magnetic axis are aligned, different 
to's become coupled and mode parity is lost. The symme- 
try with respect to S is lost even when the rotation axis 
is lined up with the magnetic field. 

2.5. Numerical aspects 

Eigenmodes and eigenvalues are calculated numerically 
using two different methods. The first method, based on 
a QZ algorithm, gives all the eigenvalues (for the discre- 
tised problem) whereas the second method is an iterative 

2 In spherical coordinates, M is defined as follows: 
MV(r,6,(fi) = V r (r,ir - 0,ip)e r -V e (r,ir - 0,<p)ee + V v (r,ir - 
0,ip)e v . 

3 In spherical coordinates, S is defined as follows: 
SV{r, 9, <p) = V r (r, 9, -<p)e r + Ve(r, 6, -tp)e g - V v (r, 9, -<p)e v . 
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Arnoldi-Chebyshev algorithm which only computes a se- 
lection of eigenvalues and their corresponding eigenmodes. 

In our numerical calculations, we use a simplified ver- 
sion of Eq. (7). Since different ra's are decoupled, we do 
not have a summation over the azimuthal order. Also, be- 
cause the divergence of both perturbations is zero, the 
functions v l m and b l m can be expressed in terms of u l m and 
a^j, respectively (see App. A of Paper I). Furthermore, in 
the case of axisymmetric modes, the summation on £ will 
actually start at 1 instead of as a result of our boundary 
conditions. In the non-axisymmetric case, the summation 
on I will start at |m| as expected. Finally, the sum on the 
spherical harmonics is truncated at L. This leads to the 
following formulas: 



v = Yl U L R T + v L s T + w L T T 

e=e min 

L 

b = J2 alRT + blST + ciT?. 



(11) 



Each of the functions u m , w m ,a m and c l m is written in the 
form of a truncated Chebyshev series: 



Nr 



k=0 

Typically, we could reach a spatial resolution of L = 
350, Nr = 930 for axisymmetric toroidal modes (when 
analysing boundary layers) and L = 858, Nr = 340 for 
non-axisymmetric modes. 

We also define Chebyshev and harmonic spectral coef- 
ficients: 



Ck 



max|u^(fc)| 
max \vL,(k)\ 



max|{4(fc)| 
iriax|u^(fc)|' 



They are useful for characterising certain aspects of a 
mode's structure, as will be seen in Sect. 3.7, as well as 
the convergence of the discretization. 

3. Axisymmetric toroidal modes 

3.1. Basic properties 

In this section, we discuss the basic properties of axisym- 
metric toroidal modes. As a reminder, the fields v and b 
are both orientated in the direction e v . From the point of 
view of spherical harmonics they are made up of (u>q, Cq) 
components. 

The equations that govern toroidal modes can be sim- 
plified. It is no longer necessary to take the curl of the 
Navier-Stokes equation to remove the pressure term, be- 
cause this term vanishes owing to axisymmetry. This leads 
to the following, nearly symmetric system: 



Xb^ = {B-V)v v - 3cos0/2r% + E m A'b, 
Xv v = (B-V)6 V + 3cosfl/2r 4 ^ + EA'v v . 



(12) 



3.2. Eigenvalue spectrum 

In many ways, the eigenvalue spectrum of toroidal 
modes is similar to that of the poloidal modes presented 
in Paper I. The complex eigenvalues are grouped into 
"horizontal" branches, as can be seen in Fig. 1. These 
"eigenbranches" , indexed by n, correspond to modes with 
n nodes in the radial direction (see Fig. 2). Both perturba- 
tions have n nodes, but in different positions. Depending 
on parity, diffusion modes (lu = 0, t ^ 0) also appear but 
these modes are of little interest to asteroseismology. In 
all that follows, we will only consider oscillatory modes. 
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Figure 1. Eigenvalue spectrum for odd modes and E = _E m = 
10~ 4 . The plusses (+) correspond to oscillatory modes and the 
diamonds (o) to diffusion modes, which only appear for the 
odd parity. The letters n and q quantify the oscillatory mode 
spectrum. 



A careful look at the eigenvalues along a single branch 
reveals a regular structure. These values are remarkably 
well lined up, much better than their poloidal counter- 
parts, and they are equally spaced out. This suggests the 
following empirical law: 



OU, 



q ■ Pn 



(13) 



Note that V • v = 0, V • b = and A' = A - 1/r 2 sin 2 



where a n and (5 n (< 0) depend on the value of n, and 
q is given in Fig. 1. By looking at the slope of the 
eigenbranches, it is possible to determine the ratio be- 
tween the real and imaginary parts of (3 n . In our case, 
Re(P n ) ~ Im(/3„). 

We can then look at the corresponding eigenmodes. 
When q increases, two phenomena appear. First of all, for 
the least-damped modes, the number of nodes in the hor- 
izontal direction seems to increase. However, these nodes 
are not all conspicuous. Secondly, modes that have a 
higher damping rate gradually depart from the magnetic 
poles. Fig. 3 illustrates both of these phenomena. Care 
must be taken with highly damped modes: the spectral 
coefficients indicate that the solutions are convergent, but 
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Figure 2. Axisymmetric toroidal modes with different vertical structures, for E — Em = 10 . The left quadrant in each image 
corresponds to magnetic dissipation and the right one to magnetic energy. From left to right, and from top to bottom, the 
modes represented here correspond to (n, q) — (2, 0), (3, 0), (4, 0) and (5, 0). A logarithmic scale is used in this figure and most 
of the other figures in this paper, as it brings out more details on the modes' structure. 



the corresponding eigenvalues are highly unstable numer- 
ically due to round-off errors. This limits the number of 
modes that can be analysed safely. 

3.3. Mode structure and resonant field lines 

The structure of the modes can be understood if one 
computes the resonance frequencies of different field lines 
for the ideal case E = E m = 0. As each field line can oscil- 
late at its own frequency, the resonance frequencies form 
a continuous spectrum called the Alfven continuum. (For 
our configuration, there will be a continuum for each value 
of n) . These eigenvalues are purely imaginary because the 
ideal MHD operator which gives the squared eigenvalues 
A 2 is self- adjoint. 

For a first approximation of these frequencies, we can 
apply the same formula as was used by Alfven (1945), 
which is derived from the time it takes for a wave traveling 
along a field line to go from one boundary to the other. 
This corresponds to a WKB-type approximation, in which 
the curvature terms have been neglected: 



-L 



di 



Field line 



\Va\ 



1 



Wf.l. 



(14) 



where n is the number of nodes along the field line. A field 
line will then be resonant if it satisfies the relation tOf \ = lo 
where to is the frequency of the mode. Since the density is 
constant, the integral can be calculated analytically and 



is given by: 



dl 



\Va\ 



sin 
1 - 



cos 7 8 3 cos 5 8 



7 



+ cos 3 8— cos 8 



if 8 V = 8 X = 



(15) 

with [/(6>)]^ = - f(6 v ). 8 V and 8 X arc the colati- 

tudes of the field line on the inner and outer boundaries 
and satisfy the relation sin 2 6 V — -q sin 2 8\ . This formula is 
the same as Eq. (26) by Kato & Watanabe (1958), except 
for the bounds on the integral. We see that the oscillation 
rate decreases for field lines further from the magnetic 
poles, hence one can expect eigenmodes with lower eigen- 
frequencies to be further from the poles. 

Taking the variations of B into account, in order to 
obtain better results, yields the following equations, which 
still only apply to individual field lines: 



ioJi.\.b = 



lUJf.l.V 



cos 8 dv 3 cos 8 



r 3 dr 
cos 8 db 



r 3 d r 



2r 4 
3cos6> 

2r 4 



v, 



(16) 



v{v) = 0, 
6(1) = 0. 



where cos 8 = \J\ — r sin 2 W\. The boundary conditions 
are derived from a simple analytical model (see Eq. (B.8)) 
and were used to establish the n + 1/2 quantization in 
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Figure 3. Different modes from the n — 3 eigenbranch. From 
top to bottom, q takes on the values 1, 3 and 9. 




Figure 4. The kinetic energy and dissipation of eigenmode 
(n,q) = (1,8) with resonant field lines superimposed. The dif- 
fusivities are E — E m = 4.10~ 5 . The dotted line corresponds to 
Eq. (14), which is a WKB-type approximation, and the dashed 
one to Eq. (16), which takes the variations of B into account. 
As opposed to other figures in this paper, a linear intensity 
scale is used instead of a logarithmic one. This gives a thinner 
appearance to the mode's structure. 



Eq. (14). With this approach, the resonant field line coin- 
cides much better with the mode (see Fig. 4). 
3.4. Polar eigenspectrum 

We will now give a precise calculation (in the adiabatic 
case (E = E m = 0)) of the resonance frequencies corre- 
sponding to the field line along the magnetic axis. This 
analysis is motivated by the role of resonating polar field 
lines in the asymptotic limit of small diffusivities. The adi- 
abatic eigenfrequencies will be denoted by the superscript 
"0" . 

These frequencies are solutions of Eq. (16) with cos 9 = 
1. By combining the two first equations and solving for v, 
we obtain: 



dr 2 



dv 
3r— + 
dr 



15 



+ uj 2 r s 



0. 



(17) 



(We have dropped the notation Wf.i..) The solution to this 
equation is (see Abramowitz & Stegun (1972)): 



Ar 2 Ji 



LOT 



+ Br 2 J_ 



uir 



(18) 



where A and B are two constants and J v is the Bessel 
function of the first kind of order v. The corresponding 
solution for b is: 



b = iBr z Ji 



cor 



iAr 2 J_ 



The boundary conditions determine the eigenfrequencies, 
via the following relation: 



LOT] 
T 



Jl 



(l) 



J_ 



torj 



J 



We will use the notation A2 = 



(20) 

2 to mean the nth 



eigenvalue calculated with this model. Asymptotically, for 
high frequencies, these values are well approximated by 
the WKB frequencies found in Paper I: w° 
The corresponding solutions are written 62 and v 
are given by the formulas: 



2(2n+l)7T 

and 



8 V 4 



+ r 2 J_ 



4 -s l 4 



n :z .J, I ^ ] J 



(21) 



4 



It is worth noting the striking similarity between these 
expressions and Eq. (30) by Kato & Watanabc (1958), 
which corresponds to the radial variation of geomagnetic 
poloidal modes (instead of toroidal ones) calculated for a 
dipolar background in which the ^-dependence of ||-B|| has 
been neglected. 

In order to show the physical meaning of these eigen- 
values, we took an eigenspectrum and prolongated the 
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n Num. 



3.8167 
17.6147 
30.6762 
43.6358 
56.5614 
69.4727 
82.3822 
95.2896 



3.8171 
17.5996 
30.6363 
43.5575 
56.4305 
69.2767 
82.1057 
94.9228 



WKB 

6.3789 
19.1367 
31.8945 
44.6524 
57.4102 
70.1680 
82.9258 
95.6836 



Table 1. A comparison between numerical polar eigenfrequen- 
cies (Num.), adiabatic ones (u>„) and a WKB approxima- 
tion. The numerical values were obtained using a least square 
method on the eigenvalues corresponding to E = Em = 10~ 4 
and correspond to the intersection between the eigenbranches 
and the imaginary axis. 



3 
<3 




Em = E 




E m =10" 3 




T 


+ 




O 



10" 5 



ID" 4 ID" 3 

Diffusivity, E 



lines formed by the eigenbranches (see Fig.l). The inter- 
sections between the lines and the imaginary axis (r = 0) 
correspond to what could be called "numerical polar eigen- 
frequencies" . The agreement between these values and A° 
is rather good, as can be seen in Tab. 1. This leads to 
the conclusion that the polar eigenfrequencies are a good 
indicator of the eigenbranches' positions. 

3.5. Asymptotic behaviour for small diffusivities 

The next step is to look at the behaviour of eigenmodes 
and eigenvalues when the diffusivities E and E m approach 
zero. There are two important reasons for this. First of all, 
realistic astrophysical values of E and E m are very small 
(< 10~ 8 ) and still out of reach for numerical solvers. Being 
able to extrapolate the behaviour of eigenmodes gives an 
educated guess as to what they would actually be for low 
diffusivities. A second reason is that it may be possible 
this way to refine analytically the adiabatic eigenvalues 
A° defined by Eq. (20) to obtain a better approximation 
of the true eigenvalues. 

We therefore analysed the behaviour of the mode n = 
1, q = for different values of the diffusivities. Figure 5 
and Fig. 6 show how its eigenvalue, its position and its 
"thickness" vary. In order to calculate the mode's position 
and thickness, we took a profile of the magnetic energy 
|6| 2 along a meridional cut with a fixed radius 4 of 0.5 and 
a horizontal extent going from 9 = to 9 = tt/2. From 
this profile, we calculated a mean value 9 and a standard 
deviation ag . These correspond to estimates of the mode's 
position and thickness, respectively. 

From these graphs it is possible to deduce some of 
the physical phenomena that arc taking place. In the case 
where E = E m , we observe that the eigenvalue approaches 
the n = 1 polar eigenfrequency given by Eq. (20). At the 
same time, the mode seems to get thinner and closer to 



Figure 5. Behaviour of eigenvalue (n, q) = (1,0) for different 
values of E and E m . Alo corresponds to the difference between 
the theoretical polar eigenfrequency, solution to Eq. (20) , and 
the actual frequency of the mode. When E m = 10~ 3 , the eigen- 
value never attains A°. When E = _E m , the slope of the Alo 
line is 1/2, meaning that Alo is proportional to E 1 ^ 2 . 
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+ 
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4 The choice of this radius is arbitrary. If we pick differ- 
ent radii, we get approximately the same empirical laws (see 
Eq. (22)). 



ID" 4 lO- 3 

Diffusivity, E 



Figure 6. Behaviour of an eigenmode for different values of E 
and E m . 8 is an estimate of the mode's position and a$ an 
estimate of the thickness. In the case E = _B m the position and 
thickness seem to be approaching zero; the slope of both lines 
is 1/4. 



the magnetic poles. We observe from Fig. 5 and Fig. 6 
that: 

Ai, - A? + £ 1 /2( r i +ioj 1 )+Et\ 

9 cx E 1 ^, (22) 

a e oc E 1 / 4 , 

in which A° is the polar eigenvalue, and r , lo 1 and r 2 are 
constants. It turns out that t 1 = lo 1 . These laws suggest 
that toroidal modes become singular when E = E m ap- 
proach zero thus revealing their true nature, namely that 
of magnetic shear layers. 
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By identifying Eq. (13) with Eq. (22), we find that: 

ai = Xl + E 1 ' 2 ^ 1 +ilo 1 )+Et 2 . (23) 
A comparison between f3\ and t 1 + iu 1 permits the iden- 



tification: 



ft =E 1 ' 2 {t 1 



(24) 



Therefore, we can cast Eq. (13) into a new form, which sep- 
arates the adiabatic contribution from the non-adiabatic 



ones: 



X n , q = A° + E l / 2 (q + \){r l n + + Et 2 



(25) 



The third term, t 2 , is roughly proportional to the square of 
the frequency and can therefore become noticeable when 
n is sufficiently large. If the eigenvalues were reduced to 
the first two terms, then for arbitrary values of E, they 
would all fall on the lines given by: 



Im(A) - uj° n _ Re (A) 



(26) 



When E m is fixed, the situation is different. The eigen- 
value does not appear to converge toward its correspond- 
ing polar eigenvalue. At the same time, the mode struc- 
ture does not become singular, which is consistent with 
the behaviour of the eigenvalue. The same experiment 
can be repeated, but with E fixed instead of E m . The 
same phenomenon appears. However, the limit eigenvalue 
for (E m = 0, E = 10~ 3 ) is not the same as for (E m = 
10- 3 ,E = 0). 

3. 6. Asymptotic solutions 

Having separated the adiabatic and non-adiabatic contri- 
butions to the eigenvalue, we will now give a justification 
of Eq. (25), based on an analytical model. We assume that 
the diffusivitics are of the form E = Ke and E m ~ K m e, 
with e approaching zero. We find that to zeroth order in e, 
eigenmodes can be put in the following form (see App. A): 

b°(r,i>) = 6°(r)/(i>), 
«°(r,f>) - w°(r)/(*), 

where (&°,v°) are gi ven by Eq. (21), and v = 
e" 1 / 4 sinO/^/r remains constant along field lines. We prove 
in App. A that the function / satisfies the following dif- 
ferential equation (sec Eq. (A. 10)): 

"d 2 .f 



(27) 



+ 



v dv v 2 



(28) 



where A° and A* are the two first orders of A, and C\, C2 
and C3 are constants given by the following formulas: 



Jr, 

C2 = /V(|6°| 2 + K J 2 )rfr, 

Jr, 

w 

Jr, 



(29) 



K m \b°J 2 + K\v°J 2 dr. 



The function / vanishes at the pole since b and v are 
toroidal. In order to mimic the confinement of eigenmodes 
near the poles, we also impose the condition that / ap- 
proaches zero as v goes to infinity. We then look for the 
eigenfunctions / and eigenvalues A* that satisfy Eq. (28). 

The variable v is rescaled to p = ve 1 ^^ \J " 2 ^f 2 where 

A° = iu>n- This leads to the following equation: 



dV + Id/ 
dp 2 p dp 



P 2+P + v /c? l C 2 C 3 /2 / 



/ = 0, (30) 



A general solution / of Eq. (30) can be expressed via 
a linear combination of the following functions: 

Sl (p) = pe-r 2 / 2 M(a,2,p 2 ), 
s 2 (p) = pe-P 2 / 2 U(a,2,p 2 ), 

a =1+ A " L % 

where M and U arc confluent hypergeometric func- 
tions, solutions to Rummer's equation (see Abramowitz 
& Stegun (1972), M is known as Rummer's function). 
When a is not a negative integer, it can be shown that 
si diverges when p goes to infinity and s 2 diverges when 
p approaches 0. Hence, the only way to satisfy our set of 
boundary conditions is to assume that a is a negative in- 
teger, namely a = —q. In that case, the functions M and 
U both become proportional to the generalised Laguerre 
polynomials (see Abramowitz & Stegun (1972)): 



I$\x) 



M{-q,2,x) = 



(x 1+q c- x ). (32) 



q+l [q + l)\dxi 
In summary, we obtain the following expressions for 
toroidal axisymmetric eigenvalues: 

-e^ 2 Xl q + 0(e), 
2(2n+ 1)tt 

(33) 



\ - X° 



A" - iloI ~ 



1 -T] 4 



\1 

A n,q 



-2(l + g)(l + i) 



Ci ' 



(34) 



The corresponding modes are given by: 
b = 6 + O(e 1 /2) j 

v = v + O(e 1 / 2 ), 

b°(r,6) = bl{r)pe-P 2 / 2 L ( q 1 \p 2 ), 
v°(r,6) = v^pe-f^L^ip 2 ), 
in which the variable p is given by the following expression: 

sin 9 



, , fir 1 1/4 

-1/4 p »tt/8 1 ^n L 2 
2C 3 



(35) 



The form of Eq. (33) corresponds perfectly to that of 
Eq. (25). A quantitative comparison reveals the accuracy 
of these formulas (see Tab. 2). It is particularly interesting 
to compare the profiles along meridional cuts with these 
theoretical predictions. Figure 7 shows such a comparison, 
for different values of q, a comparison which turns out to 
be excellent. Naturally, the analytical formulas are more 
accurate when both diffusivitics take on small values. 
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Figure 7. Comparison between different numerical profiles of the magnetic field and theoretical ones. These profiles are calculated 
along a meridional cut at a radius r ~ 0.5, and E — E m — 4.10 . 



n Num.^V) Num.ff 1 '^ 1 ) Theo. 






-0.0690 


-0.0685 


-0.0688 


1 


-0.1224 


-0.1237 


-0.1239 


2 


-0.1601 


-0.1608 


-0.1608 


3 


-0.1907 


-0.1910 


-0.1910 


4 


-0.2171 


-0.2173 


-0.2172 


5 


-0.2406 


-0.2406 


-0.2407 


6 


-0.2621 


-0.2624 


-0.2621 


7 


-0.2822 


-0.2819 


-0.2820 



Table 2. Comparison between numerical and theoretical first 
order eigenvalues (e 1 ^ 2 A 1 ). The two sets of numerical values 
(columns 2 and 3) are based on the numerical eigenspectrum 
given in Fig. 1 and are calculated using a least square method. 
These values should be very close to each other. The last col- 
umn is based on Eq. (33). 



3. 7. Boundary layers 

To complete our study of axisymmetric toroidal eigen- 
modes, we also present some results regarding the ex- 
istence of boundary layers. In our case, we expected 
the presence of Hartmann layers, which are very thin. 
Typically their non-dimensional thickness is given by 
B~ l y 'EE m in which B r is non dimensional (e.g. Potherat 
et al., 2002). This can be penalising for numerical calcu- 
lations as a high resolution is needed to properly resolve 
them. 

Two different methods were used to study boundary 
layers. The first approach consists in coming up with a 
highly simplified analytical model. In this model, devel- 
oped in App. B for clarity reasons, the permanent mag- 



netic field is constant and vertical, and the fluid domain 
is enclosed between two parallel planes. This model leads 
to the conclusion that there should be a Hartmann layer 
on the lower boundary and nothing on the upper plane. 
This Hartmann layer only produces a finite discontinuity 
of the normal gradient of both perturbations when the dif- 
fusivities vanish (see Eq. (B.3)). The model also justifies a 
posteriori the effective boundary conditions v{rj) = and 
6(1) = (see Eq. (B.8)) which were used in the calculation 
of polar eigenvalues (see Sect. 3.4). 

The second approach is numerical. By looking at the 
Chebyshev spectral coefficients, we can see the signature 
of Hartmann layers. Assuming the layer is described by 
an exponential variation as proposed by Potherat et al. 
(2002), it is possible to come up with a corresponding the- 
oretical prediction, which can then be compared with the 
actual Chebyshev spectral coefficients of the eigenmodes. 

Figure 8 shows such a comparison. The values used 
for B r are those at the magnetic poles, because of the 
mode's localization. In other words, the theoretical inner 
boundary has a thickness of rf y/ EE m and the outer one a 
thickness of y/EE m (see Eq. (1)). The theoretical bound- 
ary layers are multiplied by a constant so as to match the 
numerical ones. From this figure, we clearly deduce the 
presence of two boundary layers, one on the inside and 
one on the outside, unlike what was predicted by the ana- 
lytical model. This is so probably because in the simplified 
model, the two outer boundary conditions become redun- 
dant when E and E m approach zero, whereas they do not 
in the spherical setup. 

For larger diffusivities, the agreement between the nu- 
merical and theoretical results is not as good. This is not 
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Figure 8. Comparison between numerical (E = _E m = 10~ 4 ) 
and theoretical spectral coefficients which confirms the pres- 
ence of boundary layers. 



too surprising as the eigenmodes are in general further 
from the poles. At smaller diffusivities, the radial reso- 
lution rapidly becomes insufficient for good comparisons. 
The theoretical spectral coefficients for the inner boundary 
diminish very little for the first several hundred Chebyshev 
components. 

A final comment can be made about the numerical 
approach. By looking at the relative amplitude of the layer 
on the boundary, it is possible to determine what type of 
discontinuity it produces when both diffusivities approach 
zero. In our case, the amplitudes for both the inner and 
outer layers, obtained by comparing the numerical and 
theoretical spectral coefficients, were proportional to E in 
the case where E = E m , which suggests that there will 
be a discontinuity on the gradient of the perturbations 
when both diffusivities approach zero. This is the same 
behaviour as the single boundary layer in the analytical 
model. 



4. Non-axisymmetric modes 

We now focus on non-axisymmetric modes (m ^ 0). From 
the numerical point of view, they are twice as demand- 
ing as their axisymmetric counterparts for a given resolu- 
tion: the (u l m , aln) components are coupled to the (w^, c l m ) 
components, which gives birth to eigenvectors twice as 
large. 



haps not surprising as non-axisymmetric modes contain 
poloidal and toroidal components at the same time. 
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4.1. Eigenvalue spectrum 

As was the case for poloidal and toroidal eigenspectra, 
non-axisymmetric eigenvalues are located along "horizon- 
tal" branches. However, where there used to be one branch 
in the poloidal or toroidal case, there are now two branches 
next to each other for each m/0 (see Fig. 9). This is per- 



Figure9. Several non-axisymmetric eigenvalue spectra. The 
diffusivities E and E m both take on the value 10 -3 . These 
branches overlap with the axisymmetric ones. 



As diffusivities become small, the two branches start to 
differ. When compared with poloidal or toroidal branches, 
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it becomes immediately obvious that a first group of eigen- 
values resembles the toroidal branch and a second group 
matches its poloidal counterpart (see Fig. fO). 
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Figure 10. Comparison between axisymmetric poloidal, ax- 
isymmetric toroidal and non- axisymmetric (m = 1) eigenval- 
ues. The diffusivities E and E m both take on the values 4.10 -5 . 
The lower non-axisymmetric branch corresponds to poloidal- 
like modes. The upper branch contains toroidal-like modes 
when t is sufficiently negative. The letters "a","b","c" and "d" 
label non-axisymmetric modes that are represented in Fig. 12. 

It is then interesting to look at the structure of the 
corresponding modes. Not surprisingly, modes along the 
lower branch have a very similar appearance to that of 
poloidal modes. The radial nodes are in the same positions 
and the horizontal nodes look the same. As for the modes 
on the upper branches, it appears that the least-damped 
ones are poloidal-like, and the most damped toroidal- 
like. Figure 1 1 shows a comparison between axisymmetric 
poloidal and toroidal modes and non-axisymmetric modes 
with similar appearances. Figure 12 shows the transition 
from poloidal properties to toroidal ones, when looking at 
successive modes on the upper branch. 

One way to characterise whether a mode is more 
poloidal or toroidal is to calculate the following ratio: 



X = 



H{v! m )+H{vl)+H{al) + H{bl) 



(36) 



where H(x) represents the energy of component x. x is the 
ratio of "poloidal" energy (both kinetic and magnetic) to 
"toroidal" energy. Explicitly, it is given by: 



L 

E 

1-0 . 
■ — *-min 



vfn components whereas "toroidal" ones have strong w„ 
and c l m components (see Tab. 3). 



X |2 



l)(\vL\ 2 + \K n \ 2 )]r 2 dr 



£ /^(*+l)(|t«&| a + |4rV* 
e=e min Jr i 



(37) 

When this test is applied to the modes of Fig. 12, we 
confirm that "poloidal" modes have strong a£,, and 



Mode 


X 


a 


17.4 


b 


3.16 


c 


0.209 


d 


0.0153 



Table 3. Energy ratios as defined in Eq. (36) for modes shown 
in Fig. 10 and 12. 



4.2. Asymptotic numerical behaviour 

We can take a more detailed look at the behaviour 
of an entire pair of branches for different values of E 
and E m (keeping E — E m ), to see whether the eigen- 
values have poloidal or toroidal characteristics. Non- 
axisymmetric "poloidal" eigenvalues appear to obey the 
following empirical law: 

X = iiu" +iuj 1 E" 1 + t 1 E s , (38) 

where 7 ~ 1.5 and 5 ~ 1. This is different from what 
was observed in the case of axisymmetric toroidal modes, 
which show a slower convergence toward an asymp- 
totic adiabatic value. In addition each non-axisymmetric 
poloidal eigenvalue seems to have its own asymptotic fre- 
quency, as may also be the case for axisymmetric poloidal 
modes. A plot of the behaviour of different eigenvalues 
with respect to diffusivities is given in Fig. 13. u>» was esti- 
mated using a least-square method for the poloidal modes 
and differs for each one of them. 

Non-axisymmetric "toroidal" eigenvalues behave 
roughly like their axisymmetric counterparts. All the 
values on a single branch seem to converge toward one fre- 
quency, the corresponding polar eigenvalue (see Eq. (20)). 
This value is used in Fig. 13 as the asymptotic frequency 
of the toroidal modes. The rate at which eigenvalues 
converge is about the same as in the axisymmetric case, 
namely A, - A cx E 1 / 2 . 

It is also interesting to note that the shape of the 
toroidal eigenbranch remains more or less the same in the 
axisymmetric and non-axisymmetric cases (see Fig. 14, 
bottom). This indicates that the quantization of these 
modes in the horizontal direction is very similar to the 
one of the axisymmetric modes. However, the spacing be- 
tween each mode is not the same (see Fig. 10). 

By looking at the behaviour of the structure of differ- 
ent modes, it is possible to observe different characteristics 
(see Fig. 15). Poloidal-like modes seem to converge toward 
a non-degenerate structure. In other words, their widths 
(oe) do not approach zero, nor do their positions (0). Their 
nodes become more distinct at low diffusivities. 

Toroidal-like modes, on the other hand, exhibit more 
complex behaviour as the toroidal character appears only 
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Figure 11. Comparison between axisymmetric poloidal and toroidal modes and non-axisymmetric modes with similar properties. 
The axisymmetric modes are above, and the "corresponding" non-axisymmetric ones below. 
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Figure 13. Behaviour of r and Au; (= — w) for four 
non-axisymmetric modes. The toroidal modes start off with 
a poloidal structure, and undergo a metamorphosis into a 
toroidal appearance. The poloidal non-axisymmetric modes ap- 
proach their inviscid values much faster. 



below some given diffusivity. We understand this be- 
haviour as the consequence of the i? 1 / 4 -scale of toroidal 
modes. It is only when E 1 / 4 <C 1 that the toroidal char- 
acter may appear but this condition is not easy to satisfy 
with non-axisymmetric modes for which the lowest diffu- 
sivities reached yet are <~ 10~ 5 because of memory storage 
requirements. 



5. Discussion 

The main result of this study is that there are two basic 
behaviours for shear Alfven modes, both in the axisym- 
metric and the non-axisymmetric cases. The first one is 
the poloidal behaviour, for which modes remain regular as 
diffusivities approach zero, except for internal shear lay- 
ers located on resonant field lines. The second one is the 
toroidal behaviour: these modes become singular as dif- 
fusivities approach zero and each branch appears to have 
only one eigenfrcquency in the adiabatic limit. 

The toroidal behaviour matches the description of the 
resistive Alfven spectrum given by Kerner et al. (1986) in 
the context of tokamak-like configurations. The eigenval- 
ues lie on well-defined curves in the complex plane which 
join the endpoints of the ideal Alfven continua, as is sug- 
gested by Eq. (26). The slope of the eigenbranches is 45° 
and the damping rates of the eigenmodes are proportional 
to e 1 / 2 , which is in complete agreement with a quadratic 
variation of the permanent magnetic field ^-profile near 
the magnetic axis. 

Poloidal modes on the other hand, look like true eigen- 
modes of the ideal MHD operator with a finite number of 
resonances in the Alfven continuum (resonating field lines, 
see paper I). Such features are similar to those of "quasi- 
modes" (see Poedts & Kerner, 1991, Goossens et al., 2002) 
which combine a regular and singular behaviour; however, 
unlike quasimodes, the damping rate of our modes always 
seem to vanish in the ideal limit (e.g. Fig. 13 or Fig. 5 of 
Paper I). 
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Figure 12. Transition from poloidal to toroidal characteristics in non-axisymmetric modes. The diffusivities are E = 
E m = 4.10 -5 and m is equal to 1. The letters underneath each diagram correspond to labels on Fig. 10. 

1 /2 

weak, 0(E„i ), pressure term and are thus similar to our 
toroidal modes. 

6. Conclusions 

In this paper we completed the study of shear Alfvcn 
modes of a spherical shell started in Paper I. The global 
picture that emerges is that both non-axisymmetric and 
axisymmctric modes can be divided into two subclasses: 
the toroidal and poloidal ones. In the limit of small dif- 
fusivities, relevant for astrophysical objects, only poloidal 
modes survive since toroidal ones become singular as mag- 
netic shear layers. 

However, toroidal modes may play some important role 
in the excitation process, especially if rotation is present. 
Indeed, it is well-known (e.g. Rieutord, 1991) that the 
Coriolis force induces a coupling between toroidal and 
poloidal components, even in the axisymmetric case. This 
coupling may not strongly influence the frequency of the 
modes (except for some rotational splitting) but may rep- 
resent a channel by which modes lose energy through dis- 
sipative shear layers; in such a case this mechanism may 
decide whether a mode is or is not excited by, say, a k- 
mcchanism. 

For astrophysical applications where the compressibil- 
ity is an important aspect, pure Alfvcn modes coexist 
with magneto- acoustic modes. The incorporation of these 
modes will undoubtedly lead to a much more complicated 
spectrum of eigenvalues as well as a whole new range of 



Presently, we understand the origin of the difference 
between the toroidal and poloidal behaviour as the pres- 
ence or absence of a coupling between magnetic field lines 
by total pressure. Indeed, neglecting diffusion terms, mag- 
netic and velocity perturbations verify: 



— = - Vp* + (B • V)6 + {b ■ V)B 

8h ^ 
°-g = (B ■ V)v - (v ■ V)B 

where = p + B ■ b is the perturbed total pressure. 
These equations show that the time evolution of a per- 
turbation depends on the global form of the perturbation 
only through the pressure gradient term. This is clear if 
one considers a local perturbation, namely an Alfvcn wave 
packet. Terms like (B ■ V)6 or (B ■ V)u modulate the evo- 
lution of the perturbation according to its derivative along 
the field lines while terms like (b ■ V)J3 and (v ■ V)B in- 
troduce a dependence with respect to the local value of 
the perturbation. Last, the term Vp* makes the evolution 
of the perturbation sensitive to the pressure perturbation 
which cannot be restricted to a single field line. In pure 
toroidal modes (axisymmetric ones) this term is zero and 
therefore no coupling is possible in the ideal MHD case. In 
non-axisymmetric modes of toroidal type it is weak and 
strengthens as one shifts to modes with a strong poloidal 
character. Modes studied by Kerner et al. (1986) have a 
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Figure 15. Behaviour of non-axisymmctric modes for different 
values of the diffusivities. The "poloidal" modes seem to con- 
verge toward a non-singular structure. The "toroidal" ones, on 
the other hand, do not stop shrinking. The behaviour of their 
mean position (not shown), however, is much more irregular. 

hardly be anticipated and will be investigated in subse- 
quent work. 
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Figure 14. Behaviour of non-axisymmetric eigenbranches for 
different values of the diffusivities. E(= i? m ) takes on the val- 
ues lCT 3 ,4.1(r 4 ,2.1(n 4 ,l(r 4 and 4.KT 5 . As diffusivities di- 
minish, both eigenbranches get closer to the imaginary axis 
since |t| decreases. The first figure corresponds to "poloidal" 
eigenvalues and the second one to "toroidal" ones. The classifi- 
cation between the two types of values was done by hand. The 
dotted line represents the line given by Eq. (26) for n — 1. 



physical phenomena. In particular, we can expect new 
resonances, namely those coming from the interaction of 
eigenmodes and the slow acoustic, or cusp, continuum 
(see Sakurai et al., 1991). Taking into account compress- 
ibility will however not affect the splitting between pure 
poloidal and pure toroidal axisymmetric modes; because 
of the divergence-less nature of purely toroidal modes, we 
can expect that these modes will not be affected by com- 
pressibility and remain unchanged. All the complications 
will concentrate on the poloidal set of eigenmodes. In the 
case of non-axisymmctric modes where the poloidal and 
toroidal characters are mixed, we may expect that only 
modes with a strong poloidal shape will suffer important 
changes. However, the importance of these changes, which 
incorporate the resonances with the cusp continuum, as 
well as its influence on the shape of the spectrum, can 
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Appendix A: ^-dependence of asymptotic 
axisymmetric toroidal modes 

We begin with the toroidal equations (12), and express 
them explicitly in "semi-dipolar" coordinates (r, v, ip) . The 
coordinate v is equal to sin^/y 7 ?, and remains constant 
along field lines. Choosing v as a coordinate instead of 9 
simplifies the expression of terms like (B ■ V)b. We obtain 
the following system: 



Xb = \pY — rv 2 



Xv = vT 



rv* 



1 dv 3v 

r 3 dr 2r 4 

1 db 36 

r 3 dr 2r 4 



E m A'b, 



+ EA'v, 



(A.l) 



where A' takes on the following expression : 



(Pb 
dr 2 



' 8 2 b 



r dvdr 



iv 2 d 2 b 

4r 2 dv 2 



2db 1 

r or r a 



d 2 b 18b 
dv 2 v dv 



9v db 

4r 2 dv 

b ' 



In order to analyse the asymptotic limit of small dif- 
fusivitics, we suppose E and E m are of the form Ke and 
K m e respectively, where K and K m are two constants, and 
e a free parameter that approaches zero. We then rescale 
v to v = e~ 1 / 4 v, in accordance with our numerical results, 
and neglect terms smaller than 0(e 1 ' 2 ). This leads to the 
following system: 



Xb = (l - ^ 2 r0 2 ^j 



Xv- 



1 _ LV2 ri> 2 



1 dv 3v 

r 3 dr 2r 4 

1 db 3b 

r 3 dr 2r 4 



e 1 / 2 K 

e 1 ' 2 K 
+ — 3-9H, 



(A.2) 



where the differential operator 6 is given by: 

_ 8 2 b | 1 8b b 
dv 2 v dv v 2 

Following the form of Eq. (25), we develop the solutions 
in the following manner: 



A = A + £ V2 A i 5 
b = b° + eW, 
v = v° + e VV. 



(A.3) 



At zeroth order, we obtain the equations that give the 
polar eigenspectrum (see Eq. (16) for cos# = 1). We con- 
clude that the zeroth order solutions to Eq. (A.2) are of 
the form: 

b°(r,0) = b° n (r)f(0), 

v°{r,v) = v°„(r)f(0), (A.4) 

where (6°,w°) is given by Eq. (21), and A° is the corre- 
sponding eigenvalue. A" gives us the position of the differ- 
ent eigenbranches (as already seen in Sect. 3.4). Finding a 
quantization on the function / would enable us to explain 
the different eigenvalues on a given branch. This requires 
looking at the first order of Eq. (A.2) (in which we sub- 
stitute results from the zeroth order): 

3?! 1 



1 dv 1 

r 3 Q r 



AV-1^ 
r 3 dr 



2r 4 

351 
2r 4 



xywj b° n K m e[f] 

I Q 



2 

X° n r0 2 v° n f 



+ 



^(r = 1,0) = 0, 
v x (r = 77, z>) = 0. 

(A.5) 

(We now use the notation A* for A 1 ). This equation is 
of the form £0^1 = £1^0, where £0, £1 are respectively 
zeroth and first order linear operators, and "fo, \&i are 
zeroth and first order solutions. It is solvable only if its 
right-hand side is orthogonal to the kernel of the adjoint of 
the operator £ , with respect to the following dot product: 

r u*(r)x(r) + w*(r)y(r) 



u 




X 




w 




y _ 


H 



l*(r,0)|| ^' (A ' 6) 

where u* is the complex conjugate of u, and ||_B(r, 0) | = 
1/r 3 . The adjoint operator £j is deduced from £0 by using 
integration by parts: 



(A°)*u- 

(x ( >y w 



(A.7) 



1 dw 3w 

r 3 dr 2r 4 
1 du 3u 

r 3 dr 2r 4 - 

The adjoint boundary conditions are chosen so as to cancel 
out terms like w* (1)6(1) — w* (77)6(77), which come from the 
integration by parts: 

u(r = 1) = 0, w(r = 77) = 0. (A.8) 

The zeroth order eigenvalue A° is purely imaginary. As a 
result the operator £0 is antisymmetric: C\ = —Co. The 
boundary conditions being the same for both operators, 

0: 



the kernel of £q is the same as that of C 



u = b° n , w = v° n . (A.9) 

By calculating the dot product between (6°,i^) and 
the terms of Eq. (A.5) we obtain the following equation: 

-A^Cr/ - ^4^f + C 3 9[/] = (A.10) 
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where: 



where 5 : 



Ci 



c 2 = 



6° 

„o 



[ r 3 (\b° n \ 2 + \v° n \ 2 ) dr, 
6° 



AW + ffl dr, 

Jn 



(A.11) 



u 


1 


' K m b° n ' 


w 







A m |6' 



|2 



Equation (A. 10) gives the quantization along each eigen- 
branch. 



Appendix B: Simplified analytical model 

The model presented here justifies the use of effective 
boundary conditions throughout the paper for the adi- 
abatic calculations. In this model, we use Cartesian coor- 
dinates (x, y, z), and the fluid domain is enclosed between 
two infinite horizontal planes located at z = and z = 1. 
The permanent magnetic field is constant and in the ver- 
tical direction (B = Be z ). The perturbations take on the 
form b = b(z) e y and v = v(z)e y , in order to mimic the 
toroidal direction e v . Equation (3) becomes: 



Xb 
Xv 



dv d 2 b 

dz m dz 2 ' 

db £rd 2 i> 

dz dz 2 



(B.l) 



The corresponding boundary conditions are: 



J(o) = o, £ m £(o) = -v(0), 
J(1) = 0, 6(1) = 0. 



(B.2) 



r 



1 + (y/E - s/E^j A, 
V EE m — E, 



(B.4) 



2y/EE m 

*(c-c+) 



VEE~n 

-iX. 



The quantities fei, fc 2 , — fci,and —k 2 are solutions of the 
dispersion relation of Alfven waves in a uniform field 
(Chandrasekhar, 1961): 

-(X + Ek 2 )(X + E m k 2 ) =k 2 (B.5) 

where the non-dimensional Alfven velocity is 1 and the 
angle between the wave vector k and the permanent mag- 
netic field is 0. 

The eigenvalue A must satisfy the following relation: 



{-C~A + i[((~) 2 +^A]/cx} sinfci cos k 2 
= {-CX + i[(C~) 2 +£A]fc 2 }sinfc 2 cosfci. 



(B.6) 



If the ratio between E and E m remains constant, we 
obtain the following formula for A: 



A — — 



E + E n 



0(E 2 ), 



(B.l) 



where il n = (n + |) tt. 

When E and E m approach zero, the solutions given 
in Eq. (B.3) obey the following effective boundary condi- 
tions: 

u(0)=0 6(1) =0, (B.8) 

and their gradients admit a discontinuity on the lower 
boundary. It is interesting to note that the stress-free con- 
dition has disappeared on the lower boundary, as was ex- 
pected. On the upper boundary, it simply became redun- 
dant with the condition 6(1) = 0. 

When comparing the eigenvalues given in Eq. (B.7) to 
the ones given in Eq. (33), we notice that the term 0(e x l 2 ) 
is lacking in the first formula. This is because in the sim- 
plified model, eigenmodes are invariant in the e x direction, 
thus removing any quantization in that direction. 



Since the two conditions on the derivative of v are stress 
free conditions, they lack physical meaning when both dif- 
fusivities vanish. Therefore, we can expect them to disap- 
pear in that situation. The solutions to Eq. (B.l) are: 



v = k- 1 k 2 smk 2 (e tklZ +e lkl< - 2 -^) 

-*(l-c 2 ^)cos(fc 2 (z-l)), 

6 — k 2 sinfc 2 fc x + i£) (e 1 — e 1 ^ 5 ^ nen taking the square root of a complex number z, we 

+ (l-c 2lkl ) ((- + ^fc 2 )sin(fc 2 (z- 1)), choose the root such that arg(Vi) e] - tt/2, tt/2] 



